﻿Imports DotSpatial.Topology, DotSpatial.Data, DotSpatial.Symbology, System.Threading.Tasks
#Const SingleThreaded = True

''' <summary>
''' Encapsulates routines to convert polygon or point shapefile to raster, determine IDW weights, and reassign grid values based on attribute values
''' </summary>
''' <remarks></remarks>
Public Class clsRaster
    Private lyrFeature As FeatureLayer = Nothing
    Private FieldName As String = ""
    Private RasterName As String = ""
    Private CellSize As Double = 50
    Private NumClosestPoints As Integer = 5
    Private MaxDistance As Double = 500
    Private IDWExponent As Integer = 1
    Private lstNeighbors(,) As List(Of Tuple(Of Integer, Double)) 'indices of nearest points in feature layer and their weight (Nothing if a NoDataValue cell; 1 to NumClosestPoints otherwise)
    Private raster As Raster = Nothing

    Sub New(lyrFeature As FeatureLayer, FieldName As String, RasterFileName As String, CellSize As Double, NumClosestPoints As Integer, MaxDistance As Double, IDWExponent As Integer)
        Me.lyrFeature = lyrFeature
        Me.FieldName = FieldName
        Me.RasterName = RasterFileName
        Me.CellSize = CellSize
        Me.NumClosestPoints = NumClosestPoints
        Me.MaxDistance = MaxDistance
        Me.IDWExponent = IDWExponent
    End Sub

    Public Function CreateRaster() As Raster
        Try

            If TypeOf lyrFeature Is DotSpatial.Controls.MapPointLayer Then
                Dim rasterOptions(0) As String
                Me.raster = DotSpatial.Data.Raster.CreateRaster(RasterName, Nothing, Int(lyrFeature.DataSet.Extent.Height / CellSize) + 1, Int(lyrFeature.DataSet.Extent.Width / CellSize) + 1, 1, GetType(Double), rasterOptions)
                Me.raster.Bounds = New RasterBounds(Me.raster.NumRows, Me.raster.NumColumns, lyrFeature.DataSet.Extent)
            Else
                Me.raster = DotSpatial.Analysis.VectorToRaster.ToRaster(lyrFeature.DataSet, CellSize, FieldName, RasterName)
            End If
            Me.raster.Projection = lyrFeature.Projection
            ComputeWeights()
            UpdateGrid()
            Me.raster.Projection = ProjectProjection
            Me.raster.Save()

        Catch ex As Exception

            ErrorMsg(, ex)

        End Try

        Return Me.raster
    End Function

    Public Enum enumWContourInterpolation
        IDW
        Cressman
    End Enum

    Public Function CreateRasterWContour(WContourInterpolation As enumWContourInterpolation) As Raster
        Return Nothing
        ''note: this requires reference to WContour.dll (http://www.meteothinker.com); it is poorly documented (see chm file) so must examine demo application to see how to use it

        'Me.raster = DotSpatial.Analysis.VectorToRaster.ToRaster(lyrFeature.DataSet, CellSize, FieldName, RasterName)

        ''save outside border of this simple raster to define outline of featurelayer
        'Dim S1(1, 1) As Integer
        'Dim borders = wContour.Contour.TracingBorders(Z, X, Y, S1, .NoDataValue)

        ''discrete data array is [2,NumPoints-1] where first dimension refers to x, y, and z values
        'Dim discretepts(2, lyrFeature.DataSet.NumRows - 1) As Double
        'For i = 0 To lyrFeature.DataSet.NumRows - 1
        '    With lyrFeature.DataSet.DataTable.Rows(i)
        '        discretepts(0, i) = lyrFeature.DataSet.Features(i).Centroid.Coordinates(0).X
        '        discretepts(1, i) = lyrFeature.DataSet.Features(i).Centroid.Coordinates(0).Y
        '        discretepts(2, i) = TestNull(.Item(FieldName), Me.raster.NoDataValue)
        '    End With
        'Next
        'With lyrFeature.DataSet
        '    Dim NumRows As Integer = .Extent.Height / Me.raster.CellHeight
        '    Dim NumCols As Integer = .Extent.Width / Me.raster.CellWidth
        '    Dim X(NumCols - 1), Y(NumRows - 1) As Double
        '    For r = 0 To NumRows - 1
        '        Y(r) = .Extent.MinY + r * Me.CellSize
        '    Next
        '    For c = 0 To NumCols - 1
        '        X(c) = .Extent.MinX + c * Me.CellSize
        '    Next
        '    Dim Z(,) As Double = Nothing
        '    Select Case WContourInterpolation
        '        Case enumWContourInterpolation.IDW
        '            Z = wContour.Interpolate.Interpolation_IDW_Neighbor(discretepts, X, Y, Me.NumClosestPoints, raster.NoDataValue)
        '        Case enumWContourInterpolation.Cressman
        '            'determine average distance between closest points (polygon spacing)
        '            Dim sumDist As Double = 0.0
        '            For i = 0 To discretepts.GetUpperBound(1) - 1
        '                Dim minDist As Double = Double.MaxValue
        '                For j = 0 To discretepts.GetUpperBound(1) - 1
        '                    If i <> j Then minDist = Math.Min(mindist, Math.Sqrt((discretepts(0, i) - discretepts(0, j)) ^ 2 + (discretepts(1, i) - discretepts(1, j)) ^ 2))
        '                Next
        '                sumDist += minDist
        '            Next
        '            Dim avgDist As Double = sumDist / discretepts.GetUpperBound(1)
        '            Dim radList As New List(Of Double)
        '            For i = 10 To 2 Step -2
        '                radList.Add(i * 5) 'avgDist * 0.1)
        '            Next
        '            Z = wContour.Interpolate.Cressman(discretepts, X, Y, raster.NoDataValue, radList)
        '    End Select

        '    For r = 0 To Z.GetUpperBound(0) - 1
        '        For c = 0 To Z.GetUpperBound(1) - 1
        '            Me.raster.Value(r, c) = Z(Z.GetUpperBound(0) - r, c)
        '        Next
        '    Next
        '    Me.raster.Projection = .Projection

        '    'now trim raster using outline of feature layer
        '    Dim S1(1, 1) As Integer
        '    Dim borders = wContour.Contour.TracingBorders(Z, X, Y, S1, Me.raster.NoDataValue)

        'End With
        'raster.Save()
        'Return Me.raster
    End Function

    Public Function CreateContours(raster As Raster) As FeatureSet

        'note: this requires reference to WContour.dll (http://www.meteothinker.com); it is poorly documented (see chm file) so must examine demo application to see how to use it

        Dim Z(,), X(), Y(), Cont() As Double

        With raster
            WriteProgress(GISForm, "Computing contours...")
            ReDim Z(.NumRows, .NumColumns), X(.NumColumns), Y(.NumRows)
            For r = 0 To .NumRows - 1
                'Y(r) = (.Yllcenter + .NumRows * .CellHeight) - r * .CellHeight
                Y(r) = .Yllcenter + r * .CellHeight
            Next
            For c = 0 To raster.NumColumns - 1
                X(c) = .Xllcenter + c * .CellWidth
            Next
            For r = 0 To .NumRows - 1
                For c = 0 To raster.NumColumns - 1
                    'Z(.NumRows - r, c) = .Value(r, c)
                    Z(r, c) = .Value(r, c)
                Next
            Next
            ReDim Cont(CInt(.Maximum) - CInt(.Minimum + 1) + 1)
            For zz = 0 To Cont.Count - 1
                Cont(zz) = CInt(.Minimum + 1) + zz 'cont interval is 1
            Next
            Dim S1(1, 1) As Integer
            Dim borders = wContour.Contour.TracingBorders(Z, X, Y, S1, .NoDataValue)
            Dim contours = wContour.Contour.TracingContourLines(Z, X, Y, Cont.Count, Cont, .NoDataValue, borders, S1)

            Dim fsContour As New FeatureSet(DotSpatial.Topology.FeatureType.Line)
            fsContour.DataTable.Columns.Add("Value", GetType(Double))

            Dim yc As Double = .Yllcenter + .NumRows * .CellHeight / 2.0 'need to flip about midpoint of grid

            'For Each brd In borders
            '    For Each l In brd.LineList
            '        Dim pts As New List(Of DotSpatial.Topology.Point)
            '        For Each p In l.pointList
            '            pts.Add(New DotSpatial.Topology.Point(p.X, 2 * yc - p.Y))
            '        Next
            '        pts.Add(New DotSpatial.Topology.Point(l.pointList(0).X, 2 * yc - l.pointList(0).Y))
            '        Dim pl As New DotSpatial.Topology.LineString(pts)
            '        fsContour.AddFeature(pl)
            '    Next
            'Next
            For Each con In contours
                Dim pts As New List(Of DotSpatial.Topology.Point)
                For Each p In con.PointList
                    pts.Add(New DotSpatial.Topology.Point(p.X, 2 * yc - p.Y))
                Next
                Dim pl As New DotSpatial.Topology.LineString(pts)

                fsContour.AddFeature(pl)
                With fsContour.DataTable.Rows(fsContour.DataTable.Rows.Count - 1)
                    .Item("Value") = con.Value
                End With
            Next
            fsContour.InvalidateEnvelope()
            fsContour.InitializeVertices()
            fsContour.InitializeVertices()
            fsContour.UpdateExtent()
            fsContour.Projection = raster.Projection
            WriteProgress(GISForm)
            Return fsContour
        End With
    End Function

    Private Sub ComputeWeights()

        If lyrFeature.DataSet.NumRows = 0 Then Exit Sub

        With raster

            ReDim lstNeighbors(.NumColumns - 1, .NumRows - 1)
            Dim cnt As Integer = 0

            Dim MaxDist As Double = 0.0 'set max dist based on largest polygon

            Dim lstPoly As New List(Of Coordinate)
            For f = 0 To lyrFeature.DataSet.Features.Count - 1
                With lyrFeature.DataSet.Features(f)
                    If TypeOf lyrFeature Is DotSpatial.Controls.MapPointLayer Then
                        lstPoly.Add(.Coordinates(0))
                    Else
                        lstPoly.Add(.Centroid.Coordinates(0))
                        MaxDist = Math.Max(MaxDistance, Math.Max(.Envelope.Height, .Envelope.Width))
                    End If
                End With
            Next

            If MaxDist = 0.0 Then MaxDist = Me.MaxDistance

            Dim sw As New Stopwatch
            sw.Start()

#If SingleThreaded Then
            For r As Integer = 0 To .NumRows - 1
                WriteProgress(GISForm, cnt, .NumRows - 1, "Creating raster...")
#Else
            WriteProgress(GISForm, "Creating raster...")
            Parallel.For(0, .NumRows, Sub(r As Integer)
#End If
                                          cnt += 1

                                          For c As Integer = 0 To .NumColumns - 1
                                              If .Value(r, c) <> .NoDataValue Then
                                                  Dim x As Double = .Xllcenter + c * .CellWidth
                                                  Dim y As Double = .Yllcenter + (.NumRows - 1 - r) * .CellHeight
                                                  Dim pt As New Coordinate(x, y)
                                                  Dim dictDist As New SortedDictionary(Of Double, Integer)
                                                  For f = 0 To lstPoly.Count - 1
                                                      If dictDist.Count < NumClosestPoints Then
                                                          dictDist.Add(pt.Distance(lstPoly(f)), f)
                                                      Else
                                                          If Math.Abs(lstPoly(f).X - x) <= MaxDist AndAlso Math.Abs(lstPoly(f).Y - y) <= MaxDist Then
                                                              Dim dist As Double = pt.Distance(lstPoly(f))
                                                              If dist < MaxDist Then dictDist.Add(pt.Distance(lstPoly(f)), f)
                                                          End If
                                                      End If
                                                  Next

                                                  lstNeighbors(c, r) = New List(Of Tuple(Of Integer, Double))

                                                  With lstNeighbors(c, r)

                                                      Dim sumwt As Double = 0.0

                                                      For Each kv In dictDist
                                                          If .Count < NumClosestPoints AndAlso kv.Key < MaxDist Then
                                                              Dim wt As Double = 1.0 / kv.Key ^ Me.IDWExponent  'inverse distance raised to a power
                                                              sumwt += wt
                                                              .Add(Tuple.Create(kv.Value, wt)) 'add index and weight if closest one or within radius of influence
                                                              If .Count = NumClosestPoints Then Exit For
                                                          End If
                                                      Next

                                                      'recalculate the weights
                                                      If sumwt > 0 Then
                                                          For i = 0 To .Count - 1
                                                              .Item(i) = tuple.Create(.Item(i).Item1, .Item(i).Item2 / sumwt)
                                                          Next
                                                      End If
                                                  End With
                                              End If
                                          Next

#If SingleThreaded Then
            Next r
#Else
                                      End Sub)
#End If
            WriteProgress(GISForm)

            Debug.Print("ComputeWeights: " & sw.Elapsed.TotalSeconds)

        End With
    End Sub

    Public Sub UpdateGrid()

        If lyrFeature.DataSet.NumRows = 0 Then Exit Sub

        With raster

            Dim sw As New Stopwatch
            sw.Start()

#If SingleThreaded Then
            For r As Integer = 0 To .NumRows - 1
#Else
            Parallel.For(0, .NumRows, Sub(r As Integer)
#End If

                                          For c As Integer = 0 To .NumColumns - 1
                                              Dim v As Double = 0
                                              If lstNeighbors(c, r) Is Nothing Then
                                                  v = .NoDataValue
                                              Else
                                                  For Each t In lstNeighbors(c, r) 'item 1 is index, 2 is weight
                                                      Dim f As Feature = lyrFeature.DataSet.Features(t.Item1)
                                                      v += TestNull(f.DataRow.Item(FieldName), 0) * t.Item2
                                                  Next
                                              End If
                                              .Value(r, c) = v
                                          Next

#If SingleThreaded Then
            Next r
#Else
                                      End Sub)
#End If

            Debug.Print("UpdateGrid: " & sw.Elapsed.TotalSeconds)
        End With
    End Sub

    ''' <summary>
    ''' Create coloring scheme for raster layer based on current values
    ''' </summary>
    Friend Shared Sub CreateScheme(lyrRaster As RasterLayer)
        Try
            With GISForm.AnimationInfo

                Dim rs As New RasterSymbolizer
                Dim cs As New ColorScheme
                cs.SuspendChangeEvent()
                cs.SuspendEvents()

                If .OverrideScheme Then
                    'For i As Integer = 0 To GISForm.AnimationInfo.CustomScheme.Count - 1
                    '    With GISForm.AnimationInfo.CustomScheme(i)
                    '        Dim fillcolor As Color = .Item2
                    '        Dim pc As PolygonCategory

                    '        'create new category for symbology
                    '        If i = 0 Then
                    '            pc = New PolygonCategory(GISForm.AnimationInfo.Colors(0), Color.Gray, 0)
                    '            pc.FilterExpression = String.Format("[{0}]<={1:0.000}", "RESULT", .Item1)
                    '            pc.LegendText = String.Format("<{0:0.000}", .Item1)
                    '            polyScheme.AddCategory(pc)
                    '        End If
                    '        pc = New PolygonCategory(fillcolor, Color.Gray, 0)
                    '        If i = GISForm.AnimationInfo.CustomScheme.Count - 1 Then
                    '            pc.FilterExpression = String.Format("[{0}]>={1:0.000}", "RESULT", .Item1)
                    '            pc.LegendText = String.Format(">{0:0.000}", .Item1)
                    '        Else
                    '            pc.FilterExpression = String.Format("[{0}]>={1:0.000} and [{0}]<={2:0.000}", "RESULT", .Item1, GISForm.AnimationInfo.CustomScheme(i + 1).Item1)
                    '            pc.LegendText = String.Format("{0:0.000} - {1:0.000}", .Item1, GISForm.AnimationInfo.CustomScheme(i + 1).Item1)
                    '        End If
                    '        'add category to scheme
                    '        polyScheme.AddCategory(pc)
                    '    End With
                    'Next
                Else
                    'create three categories using 3 defined colors in options (color 0 is for undefined)

                    Dim NumEntries As Integer = 3

                    Dim minValue As Double = .MinValue
                    Dim maxValue As Double = .MaxValue
                    If minValue = maxValue AndAlso minValue = 0 Then
                        minValue = lyrRaster.Minimum
                        maxValue = lyrRaster.Maximum
                    End If

                    For i As Int16 = 0 To NumEntries - 1

                        'create new category for symbology
                        Dim range As Double = (maxValue - minValue) / NumEntries
                        Dim start As Double = minValue + range * i
                        Dim [end] As Double = start + range
                        Dim rc As New ColorCategory(start, [end], .Colors(i + 1), .Colors(i + 2))
                        rc.LegendText = String.Format("{0:N2} - {1:N2}", start, [end])

                        'add category to scheme
                        cs.AddCategory(rc)

                    Next
                End If

                cs.ResumeChangeEvent()
                cs.ResumeEvents()

                rs.Scheme = cs

                'set scheme to layer
                lyrRaster.Symbolizer = rs

                GISForm.GISMap.Refresh()

            End With

        Catch ex As Exception
            ErrorMsg(, ex)
        End Try
    End Sub

    ''' <summary>
    ''' Set all raster grid cells outside of polygon to nodatavalue
    ''' </summary>
    ''' <param name="Poly">Clipping polygon</param>
    ''' <remarks></remarks>
    Friend Sub ClipWithPolygon(Poly As Polygon)
        With Me.raster

#If SingleThreaded Then
            For r = 0 To .NumRows - 1
#Else
                    Parallel.For(0, .NumBands, Sub(r)
#End If

                For c = 0 To .NumColumns - 1
                    Dim pt As New Point(.Xllcenter + c * .CellWidth, .Yllcenter + r * .CellHeight)
                    If .Value(r, c) <> .NoDataValue AndAlso Not Poly.Contains(pt) Then
                        .Value(r, c) = .NoDataValue
                    End If
                Next

#If SingleThreaded Then
            Next
#Else
                                               End Sub)
#End If
        End With
    End Sub

    ''' <summary>
    ''' Set all raster grid cells outside of all polygons to NoDataValue
    ''' </summary>
    ''' <param name="lyrRaster">Raster layer to clip</param>
    ''' <param name="lyrPoly">Clipping polygon feature layer</param>
    ''' <remarks></remarks>
    Friend Shared Sub ClipWithPolygon(lyrRaster As RasterLayer, lyrPoly As FeatureLayer)
        Try

            With lyrRaster.DataSet

#If SingleThreaded Then
                For r = 0 To .NumRows - 1
                    WriteProgress(GISForm, r, .NumRows - 1, "Clipping raster...")
#Else
                Parallel.For(0, .NumRows, Sub(r)
#End If
                    For c = 0 To .NumColumns - 1
                        Dim pt As New Point(.Xllcenter + c * .CellWidth, .Yllcenter + (.NumRows - r - 1) * .CellHeight)
                        Dim found As Boolean = False
                        For Each f In lyrPoly.DataSet.Features
                            If .Value(r, c) <> .NoDataValue AndAlso CType(f.BasicGeometry, Polygon).Contains(pt) Then found = True : Exit For
                        Next
                        If Not found Then .Value(r, c) = .NoDataValue
                    Next
#If SingleThreaded Then
                Next
#Else
                                          End Sub)
#End If
                .Save()
            End With

            'after changing values, cannot figure out how to refresh the raster; caller will need to remove and reload the raster and reapply the color scheme

        Catch ex As Exception

            ErrorMsg(, ex)

        Finally

            WriteProgress(GISForm)

        End Try

    End Sub

End Class
